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HIGHLIGHTS 


►  New  boundary  conditions  were 
derived  to  account  for  the  Stern 
layer  without  simulating  formally  it. 

►  Both  the  Stern  and  diffuse  layers 
were  rigorously  accounted  for. 

►  The  new  boundary  conditions  were 
valid  for  planar,  cylindrical,  and 
spherical  electrodes  or  pores. 

►  They  made  possible  the  simulations 
of  EDLCs  with  highly  ordered  porous 
electrodes. 

►  Predictions  of  capacitances  of  an 
EDLC  with  ordered  bimodal  meso- 
porous  carbons  agreed  well  with 
experimental  results. 
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This  paper  presents  general  mathematical  formulations  for  simulating  electric  double  layer  capacitors 
(EDLCs)  with  three-dimensional  ordered  structures.  For  the  first  time,  a  general  set  of  boundary  condi¬ 
tions  was  derived  in  order  to  account  for  the  Stem  layer  without  simulating  it  in  the  computational 
domain.  These  boundary  conditions  were  valid  for  planar,  cylindrical,  and  spherical  electrode  particles  or 
pores.  They  made  possible  the  simulations  of  EDLCs  with  complex  geometries  while  rigorously 
accounting  for  both  the  Stern  and  diffuse  layers.  The  model  also  simultaneously  accounted  not  only  for  3D 
electrode  morphology  but  also  for  finite  ion  size  and  field-dependent  electrolyte  dielectric  permittivity.  It 
was  used  to  faithfully  simulate  the  complex  structure  of  an  EDLC  electrode  consisting  of  ordered  bimodal 
mesoporous  carbon  featuring  both  macropores  and  mesopores.  Areal  and  gravimetric  capacitances  were 
predicted  based  on  non-solvated  and  solvated  ion  diameters.  These  two  cases  set  the  upper  and  lower 
bounds  for  the  predicted  capacitances.  The  capacitances  predicted  using  non-solvated  ion  diameter  were 
found  to  be  in  good  agreement  with  experimental  data  reported  in  the  literature.  All  surfaces  contributed 
to  the  overall  capacitance  of  EDLCs.  The  gravimetric  capacitance  of  different  bimodal  carbons  increased 
linearly  with  increasing  specific  surface  area  corresponding  to  constant  areal  capacitance. 
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1.  Introduction 

Electric  double  layer  capacitors  (EDLCs)  have  been  the  subject  of 
intense  studies  in  recent  years  due  to  their  promises  as  electrical 
energy  storage  devices  [1—5],  EDLCs  store  electric  charges  physi¬ 
cally  in  the  electric  double  layer  forming  at  electrode/electrolyte 


interfaces  accessible  to  ions  present  in  the  electrolyte  [1 — 5]. 
Fig.  1(a)  shows  a  schematic  of  the  electric  double  layer  structure 
forming  near  a  positively  charged  cylindrical  or  spherical  electrode 
particle.  Solvated  anions  of  diameter  a  migrate  and  adsorb  to  the 
electrode  surface  due  to  electrostatic  forces  while  the  cations  are 
repelled  [6—9],  The  Stem  layer  is  defined  as  the  compact  layer  or 
inner  layer  near  the  electrode  surface  [6—9].  Note  that  there  are 
no  free  charges  within  the  Stern  layer  [6-8],  Beyond  the  Stern 
layer  is  the  so-called  diffuse  layer  where  anions  and  cations  are 
mobile  under  the  coupled  influence  of  electrostatic  forces  and 
diffusion  [6-9], 

Electrodes  in  EDLCs  are  typically  made  of  mesoporous  materials 
offering  large  surface  area.  Research  efforts  have  focused  on 
increasing  the  energy  and  power  densities  of  supercapacitors  by 
increasing  the  surface  area  of  porous  electrodes  and  tailoring  their 
morphology  or  pore  size  distribution  [  1  — 5  ] .  In  particular,  electrodes 
with  three-dimensional  ordered  structures  [10-27]  have  attracted 
significant  attention  due  to  [23-29]  (1)  their  small  ion  transport 


resistance,  (2)  their  uniform  pore  connection  leading  to  short  ion 
diffusion  length,  and  (3)  their  continuous  electron  transport 
framework.  For  example,  Woo  et  al.  [  1 5  ]  synthesized  highly-ordered 
mesoporous  carbon  films  as  electrodes  for  EDLCs.  These  carbon 
films  had  ordered  “bimodal”  structure  featuring  both  inter¬ 
connected  macropores  and  mesopores.  In  particular,  their  “CP204- 
S15”  carbon  film  had  specific  surface  area  of  Sbet  =  1003  m2  g  1 
[15],  The  radius  of  the  macropores  and  mesopores  was  reported  to 
be  95  nm  and  7  nm  [15],  respectively.  The  surface  area  due  to 
micropores  with  diameter  around  or  less  than  2  nm  was  less  than  6% 
[15].  From  the  FE-SEM  image  (Fig.  3(a)  in  Ref.  [15]),  the  radius  of  the 
channels  between  macropores  in  “CP204-S15”  carbon  film  was 
estimated  to  be  20  nm  while  the  carbon  wall  thickness  was  about 
2  nm.  The  electrolyte  was  1  mol  L-1  (C2H5)4NBF4  (or  TEABF4)  in 
propylene  carbonate  while  the  potential  window  was  2  V.  The 
capacitances  were  measured  using  galvanostatic  charge/discharge 
at  low  current  density  40  mA  g-1  using  the  three-electrode 
method  [15],  Finally,  the  areal  and  gravimetric  capacitances  for 
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Cylindrical  electrode  radius,  R0  (nm) 


Cylindrical  pore  radius,  R0  (nm) 


Fig.  3.  Predicted  Stern  layer,  diffuse  layer,  and  total  areal  capacitances  for  (a)  a  cylin¬ 
drical  electrode  particle  and  (b)  a  cylindrical  pore  as  a  function  of  radius  R0  ranging 
from  2.5  to  60  nm.  Results  were  obtained  using  conventional  and  our  new  boundary 
conditions  (BC)  with  a  =  0.68  nm,  c»  =  1  mol  IT1,  (,  =  2  V,  and  electrolyte 
permittivity  given  by  Equation  (5). 


the  “CP204-S15”  carbon  film  were  reported  to  be  Cs  -  9.4  pF  cm  2 
and  Cg  =  95  F  g_1,  respectively  [15]. 

Numerous  experimental  studies  have  been  devoted  to  charac¬ 
terizing  the  performances  of  EDLCs  and  assessing  the  effects  of 
electrode  morphology  as  well  as  of  the  physical  or  electrochemical 
properties  of  electrodes  and  electrolytes  [1-5,10—27],  Experimental 
approaches  are  typically  time  consuming  and  costly.  They  also  rely 
on  trial  and  errors  in  order  to  optimize  EDLCs.  On  the  contrary, 
accurate  and  reliable  numerical  tools  can  facilitate  the  design  and 
optimization  of  the  electrode  morphology  in  a  more  systematic  and 
efficient  way.  Moreover,  they  can  be  used  to  understand  the  elec¬ 
trochemical  and  transport  processes  involved  in  EDLCs  [30,31  ].  For 
example,  they  can  predict  the  local  electric  potential  and  ion 
concentrations  throughout  the  mesoporous  electrodes  [32-34] 
which  cannot  be  measured  experimentally.  However,  such  numer¬ 
ical  simulations  are  complicated  by  the  multi-scales  (from  sub¬ 
nanometer  to  micron)  and  multi-physics  nature  of  the  problem. 
They  should  also  be  validated  against  experimental  data. 

This  paper  aims  to  develop  a  three-dimensional  (3D)  model 
based  on  continuum  theory  for  simulating  EDLCs  with  ordered 
mesoporous  electrode  structures.  For  the  first  time,  the  model 
simultaneously  and  rigorously  accounts  for  (1)  3D  electrode 
morphology,  (2)  finite  ion  size,  (3)  Stern  and  diffuse  layers,  and  (4) 


the  dependency  of  the  electrolyte  dielectric  permittivity  on  the 
local  electric  field.  First,  a  new  set  of  boundary  conditions  was 
derived  to  account  for  the  Stern  layer  without  simulating  it  in  the 
electrolyte  domain.  The  model  was  then  used  to  simulate  faithfully 
the  electrode  morphology  of  CP204-S15  mesoporous  carbon  EDLC 
synthesized  and  characterized  by  Woo  et  al.  [15]. 


2.  Background 

2.1.  Traditional  modeling  approaches 


Equivalent  RC  circuit  models  and  more  complex  transmission 
line  models  [35—37]  have  been  traditionally  used  to  simulate 
EDLCs.  However,  these  models  inherently  neglect  ion  diffusion  and 
non-uniform  ion  concentration  in  the  electrolyte  [30,31,38-40]. 
Thus,  these  models  may  not  be  valid  for  EDLCs  under  large  electric 
potential  and  electrolyte  concentration  [30,31,38—41],  Alterna¬ 
tively,  homogeneous  models  have  also  been  developed  to  investi¬ 
gate  the  charging/discharging  dynamics  of  EDLCs  [42-47],  These 
models  treat  the  heterogeneous  microstructure  of  the  electrodes  as 
homogeneous  with  some  effective  macroscopic  properties  deter¬ 
mined  from  effective  medium  approximations  and  depending  on 
porosity  and  specific  area  [42—47],  In  addition,  they  typically 
impose  the  areal  capacitance  or  volumetric  capacitance  rather  than 
predict  them  [42-47], 

Moreover,  the  Helmholtz  model  [48—50]  and  Gouy— Chapman- 
Stem  (GCS)  model  [51]  are  frequently  used  to  simulate  EDLCs  with 
one-  or  two-dimensional  electrode  structure.  In  these  models,  the 
electrolyte  dielectric  permittivity  is  either  assumed  to  be  constant 
[50,51]  or  treated  as  a  fitting  parameter  [48-50].  However,  the 
relative  permittivity  er  of  polar  electrolytes  is  known  to  significantly 
decrease  as  the  electric  field  increases  [52—54],  In  addition,  the  GCS 
model  neglects  the  finite  size  of  ions  and  treat  ions  as  point-charges 
[32-34,55,56].  This  assumption  breaks  down  when  either 
the  electrolyte  concentration  c„  or  the  electric  potential  is  large 
[30—34,55,56],  Therefore,  the  GCS  model  is  invalid  for  practical 
EDLCs  with  typical  electrolyte  concentration  larger  than  1  mol  L-1 
and  potential  window  larger  than  1  V  [34], 

Due  to  their  intrinsic  limitations,  none  of  the  above-mentioned 
models  can  account  for  the  three-dimensional  mesoporous  elec¬ 
trode  morphology.  The  first  equilibrium  simulations  of  EDLCs  with 
three-dimensional  electrode  morphology  were  reported  by  Pilon 
and  co-workers  [32,33].  These  simulations  also  accounted  for  finite 
ion  size  as  well  as  the  dependency  of  the  electrolyte  dielectric 
permittivity  on  the  local  electric  field  [32,33].  However,  the 
computations  of  the  Stem  and  diffuse  layer  capacitances  were 
decoupled  due  to  the  complex  electrode  structures  [32,33].  Our 
recent  study  indicated  that  the  Stern  and  diffuse  layer  need  to  be 
simulated  simultaneously  in  order  to  predict  accurately  the  electric 
double  layer  capacitances  [34],  To  the  best  of  our  knowledge,  the 
Stem  and  diffuse  layers  have  been  simultaneously  simulated  only 
for  one-  or  two-dimensional  electrode  structures  such  as  planar 
electrodes  [30,31]  and  a  single  cylindrical  or  spherical  electrode 
particle  or  pore  [34,48-51  ]. 

Under  equilibrium  conditions,  the  local  electric  potential  f{r)  at 
location  r  in  the  electrolyte  can  be  found  by  solving  the  modified 
Poisson— Boltzmann  (MPB)  model  with  a  Stern  layer  accounting  for 
the  finite  ion  size  and  expressed  as  [34,55,56], 


V  •  (fQfrVi/')  =  0  in  the  Stern  layer 
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Then,  the  local  ion  concentration  c(r)  is  given  by  [55,56] 


c(r) 


Cooexp(  -z0e\l//kBT) 

’+2”“2(S) 


(2) 


where  eo  and  er  are  the  free  space  permittivity 
(eo  —  8.854  x  10-12  F  m-1)  and  the  relative  permittivity  of  the 
electrolyte  solution,  respectively.  The  valency  of  the  symmetric 
electrolyte  is  denoted  by  Zo,  while  c»  is  the  bulk  molar  concen¬ 
tration  of  electrolyte,  T  is  the  absolute  temperature,  e  is  the 
elementary  charge  (e  =  1.602  x  10-19  C),  N&  and  kg  are  the  Avo- 
gadro’s  number  ( Na  —  6.022  x  1023  mol-1)  and  the  Boltzmann 
constant  (kg  =  1.381  x  10-23  m2kg  K-1  s-2),  respectively.  The 
packing  parameter  is  defined  as  v  =  2a3N>\C„  where  a  is  the  effec¬ 
tive  ion  diameter.  It  represents  the  ratio  of  the  total  bulk  ion 
concentration  to  the  maximum  ion  concentration  cm  =  1  //\Ua3 
assuming  a  simple  cubic  ion  packing  [40,55,56],  Therefore,  v  should 
not  be  larger  than  unity  for  the  model  to  be  physically  acceptable. 


=  Hj  =  Cf  [*,  -  =  H)]  (4) 

where  Cf1  =  f0er/H  is  the  Stem  layer  capacitance  for  planar  elec¬ 
trodes  [38,40,55,56,59—66],  Then,  Equations  (3c)  and  (4)  form 
a  complete  set  of  boundary  conditions  for  the  entire  electric  double 
layer  while  simulating  only  the  diffuse  layer  from  x  =  H  to  x  =  L,.  To 
the  best  of  our  knowledge,  no  similar  approach  has  been  proposed 
for  simulating  electric  double  layers  near  electrodes  in  other 
geometries  or  coordinate  systems. 

This  paper  aims  to  develop  a  new  set  of  boundary  conditions  to 
account  for  the  Stern  layer  without  simulating  it  in  the  computational 
domain  for  planar,  cylindrical,  and  spherical  electrode  particles  or 
pores.  This  presents  the  advantages  of  simplifying  the  computational 
domain  by  simulating  only  the  diffuse  layer  and  thus  reducing  the 
number  of  finite  elements  and  the  computational  cost  and  time. 
Moreover,  it  also  enables  the  simulations  of  three-dimensional 
highly-ordered  mesoporous  electrode  structures.  Finally,  this 
approach  was  demonstrated  and  validated  against  experimental  data 
for  ordered  bimodal  carbon  films  reported  in  Ref.  [15]. 


2.2.  Conventional  boundary  conditions 


Boundary  conditions  are  required  to  predict  the  electric  poten¬ 
tial  and  ion  concentration  profile  in  the  electrolyte.  The  electric 
potential  at  the  electrode/electrolyte  interface  is  typically 
prescribed  under  equilibrium  conditions  [6-9,32,34,56,57],  For 
a  sphere  or  cylinder  of  radius  Ro,  it  is  given  by 

iKr  =  Ro)  =  fs,  (3a) 

In  addition,  the  electric  potential  and  displacement  are  contin¬ 
uous  across  the  Stern/diffuse  layer  interface  located  at  r  =  Ro  4-  H  so 
that  [6,30,31,58], 


t(r  =  R0  +  H =V'(r  =  Ro  +  H+)  and 
eoer^  (r  =  R0+  H-)  =  e0er^  (r  =  Ro  +  H+) 


Far  away  from  the  electrode  surface,  the  electric  potential  and 
ion  concentration  are  constant  such  that  [30—34,38], 


<A(r  =  R0  +Ij)  =0  and  c,-(r  =  R0  +L,)  =  c«,  (3c) 

In  fact,  the  presence  of  the  very  thin  Stem  layer  near  the  electrode 
surface  causes  several  numerical  challenges.  First,  the  Stern  layer 
complicates  the  computational  domain  by  introducing  an  additional 
length  scale  which  is  significantly  smaller  than  that  of  the  diffuse 
layer.  Therefore,  the  computational  domain  becomes  extremely 
complicated  and  the  number  of  meshes  prohibitively  large  for 
simulating  three-dimensional  electrode  structures.  Second,  the 
governing  equations  for  the  electric  potential  and  ion  concentrations 
in  the  Stern  and  diffuse  layers  are  numerically  solved  separately  and 
coupled  through  the  boundary  conditions  [Equation  (3b)].  These 
equations  must  be  solved  simultaneously  thus  requiring  excessive 
computational  time  and  resources.  Therefore,  the  MPB  model  with 
a  Stern  layer  [Equation  (1)]  and  the  conventional  boundary  condi¬ 
tions  [Equation  (3)]  [30,31,34,51]  cannot  be  used  to  simulate  three- 
dimensional  structures  such  as  those  encountered  in  practical 
EDLCs. 

Alternatively,  the  Stem  layer  forming  near  planar  electrodes  can 
also  be  accounted  for  via  a  modified  boundary  condition  without 
simulating  it  explicitly  in  the  electrolyte  domain.  In  one¬ 
dimensional  Cartesian  coordinates,  the  corresponding  boundary 
condition  at  the  Stern/diffuse  layer  interface  located  at  x  =  H  has 
been  derived  as  [38,40,55,56,59-66], 


3.  Analysis 

3.1.  Schematics  and  assumptions 

Fig.  2(a)  shows  the  schematic  representation  of  the  ordered 
bimodal  mesoporous  carbon  simulated  in  this  study.  Here,  the 
dimensions  of  the  simulated  electrode  structure  were  identical  to 
those  of  the  bimodal  mesoporous  carbon  “CP204-S15”  reported  in 
Ref.  [15]  as  previously  discussed.  Note  that  the  capacitances  pre¬ 
dicted  by  the  MPB  model  with  the  Stern  layer  [Equation  (1)]  are 
identical  for  the  positive  and  negative  electrodes  in  binary  and 
symmetric  electrolytes  [32,34,55,56]  such  as  those  considered  here. 
Therefore,  it  suffices  to  simulate  only  one  electrode  [32,34,55,56]. 
In  addition,  the  present  study  simulated  one  unit  cell  of  the  3D 
porous  electrode  structure.  Further  increasing  the  number  of  unit 
cells  was  found  to  have  no  effect  on  the  predicted  areal  and 
gravimetric  capacitances  under  equilibrium  conditions  [32],  By 
virtue  of  symmetry,  it  suffices  to  simulate  only  1  /12th  of  one  unit 
cell.  Fig.  2(b)  shows  the  schematic  of  the  computational  domain 
simulated  in  this  study.  The  density  of  carbon  materials  is  about 
p  =  1.8  g  cm-3  in  its  amorphous  phase  [67],  Then,  the  number  of 
mesopores  existing  in  the  walls  was  adjusted  so  that  the  specific 
area  of  the  simulated  electrode  matched  that  of  the  actual  bimodal 
electrodes  ranging  from  492  to  1504  m2  g-1  [15],  Overall,  the 
specific  area  of  the  simulated  structure  [Fig.  2(b)]  was  adjusted  to 
be  960  m2  g-1  which  falls  within  5%  of  its  experimental  value  of 
1003  m2  g-1  [15].  The  thickness  of  the  electrolyte  region  at  the  edge 
of  the  electrode  was  specified  to  be  Li  =  30  nm.  Increasing  this 
thickness  to  L,  =  60  nm  was  found  to  have  no  effect  on  the  predicted 
areal  and  gravimetric  capacitances  due  to  the  rapid  decrease  of  the 
electric  potential  from  the  electrode  surface  caused  by  the  thin 
electric  double  layer  [32—34], 

To  make  the  problem  mathematically  tractable,  the  following 
assumptions  were  made:  (1)  the  electric  potential  and  ion 
concentrations  reached  their  equilibrium  states,  (2)  anions  and 
cations  had  the  same  effective  diameter  assumed  to  be  constant  and 
independent  of  electrolyte  concentration  [55,56,68],  (3)  the  elec¬ 
trolyte  relative  permittivity  was  constant  and  uniform  within  the 
Stern  layer.  Note  that  the  electrolyte  relative  permittivity  is  typically 
defined  for  media  with  characteristic  length  larger  than  1  or  2  nm 
[69,70],  (4)  isothermal  conditions  prevailed  throughout  the  elec¬ 
trode  and  electrolyte,  (5)  advection  of  the  electrolyte  was  assumed 
to  be  negligible,  (6)  the  ions  could  only  accumulate  at  the  Stern/ 
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diffuse  layer  interface  and  could  not  diffuse  into  the  electrode,  i.e., 
there  was  no  ion  insertion,  and  (7)  the  specific  ion  adsorption  due 
to  non-electrostatic  forces  were  assumed  to  be  negligible. 

3.2.  Constitutive  relations 


Equation  (8)  relates  the  local  electric  potential  to  its  gradient  at 
r  =  Ro  +  H.  It  serves  as  a  new  boundary  condition  at  the  Stem / 
diffuse  layer  interface. 

Similarly,  for  spherical  electrodes,  the  boundary  condition  at  the 
Stem/diffuse  layer  interface  can  be  derived  as, 


In  order  to  solve  Equations  (1 )  and  (3),  the  electrolyte  properties 
er,  Zo,  c„o,  and  a  along  with  the  temperature  T  and  the  surface 
potential  \j/s  are  needed.  In  the  diffuse  layer,  the  Booth  model  was 
used  to  account  for  the  dependency  of  electrolyte  dielectric 
permittivity  on  the  local  electric  field  [32,34,52—54], 

er(E)  =n2  +  (er(0)  -  n2)  jcoth(/?E)  -  for 

E  >  107V  m_1  (5a) 


-w,(£„)|(r,Ro+H) 

xW'j-i Kr  =  Ro  +  H)}  (9) 

Note  that  Equations  (4),  (8)  and  (9)  can  be  rewritten  in  a  gener¬ 
alized  form  for  planar,  cylindrical,  and  spherical  electrodes  as, 

-coerCEHW-  (g)  =  C» (Ro%)Vs  -  Hrn)]  (10) 


er  =  er(0)  for  E<107Vm1  (5b) 

where  £  =  |—V^|  is  the  norm  of  the  local  electrical  field  vector,  er(0) 
is  the  relative  permittivity  at  zero  electric  field,  and  n  is  the  index  of 
refraction  of  the  electrolyte  at  zero  electric  field  frequency.  In  the 
Stern  layer,  the  electrolyte  dielectric  permittivity  was  imposed  as 
constant  and  uniform  [assumption  (3)].  Its  value  was  evaluated 
using  Equation  (5)  and  the  computed  local  electric  field  at  the 
Stern/diffuse  layer  interface. 

The  electrolyte  solution  used  in  Ref.  [15]  was  TEABF4  in 
propylene  carbonate  solution  at  room  temperature  (T  =  298  K) 
characterized  by  the  following  properties:  rfO)  =  64.4  [71  ],  n  —  1.42 
[72],  /?  =  1.314  x  10-8  m  V-1  [32],  and  Zo  —  1.  The  ion  diameter  of 
non-solvated  TEA+  and  BF4  ions  are  a  =  0.68  and  0.34  nm  [73], 
respectively.  Their  solvated  ion  diameters  were  reported  to  be 
a  =  1.36  and  1.40  nm  [74,75],  respectively.  The  electrolyte 
concentration  and  the  surface  electric  potential  were  chosen  to  be 
the  same  as  those  used  in  Ref.  [15],  namely,  c»  =  1.0  mol  L-1  and 
i/'s  =  2  V.  In  addition,  the  Stern  layer  thickness  H  was  approximated 
as  the  radius  of  solvated  ions,  i.e.,  H  =  a/2  [7-9], 

3.3.  Derivation  of  generalized  boundary  conditions 

This  section  presents  the  derivation  of  a  generalized  boundary 
condition  valid  for  cylindrical  and  spherical  electrode  particles  or 
pores.  For  a  cylindrical  electrode  of  radius  Ro  (Fig.  1(a)),  Equation  (1 ) 
in  the  Stern  layer  is  expressed  as  [30,31,51  ] 

3 \  =0  for  R0<r<Ro  +  H  (6) 

where  er  =  ef_En)  is  the  uniform  electrolyte  relative  permittivity 
within  the  Stern  layer  [assumption  (3)].  Its  value  is  evaluated  at  the 
Stern/diffuse  layer  interface  located  at  Th  =  Ro  +  H  using  Booth 
model  [Equation  (5)]  based  on  the  local  electric  field  E/^rw).  Then, 
integrating  Equation  (6)  twice  with  respect  to  r  from  r  =  Ro  to 
r  =  Ro  +  H  using  the  boundary  condition  given  by  Equations  (3a) 
and  (3b)  yields, 


where  rn  is  the  local  position  vector  at  the  Stern/diffuse  layer 
interface  located  at  r^  =  Ro  +  H  for  cylindrical  and  spherical  elec¬ 
trodes.  Note  that  r h/ch  represents  the  local  outward  normal  unit 
vector  at  the  Stern/diffuse  layer  interface.  Here,  Cf  is  the  Stern  layer 
capacitance  predicted  by  the  Helmholtz  model  assuming  constant 
er  within  the  Stern  layer  and  given  by  [34,76], 

_  eofr(EH)  for  planar  electrode  11a 

C?  =  (for  cylindrical  electrode  lib 

s  R0log(l  +  H/R0) 

ch  for  spherical  electrode  11c 

The  value  of  p  in  Equation  (10)  is  given  by, 
p  =  0  for  planar  electrodes  12a 

p  =  l  for  cylindrical  electrodes  12b 

p  =  2  for  spherical  electrodes  12c 

Moreover,  for  cylindrical  and  spherical  pores  of  radius  R0  illus¬ 
trated  in  Fig.  1(b),  the  new  boundary  condition  at  the  Stem/diffuse 
layer  interface  located  atru  =  R0  -  H  can  be  similarly  derived  as, 

-eoCrC EHW-  (*£j  =  C»  (^)^s  -  i(rH)]  (13) 

where  p  =  1  or  2  for  cylindrical  or  spherical  pores,  respectively. 
Here  also,  the  Stern  layer  capacitance  C''  for  cylindrical  or  spherical 
pores  is  given  by  the  Helmholtz  model  assuming  constant  er  within 
the  Stern  layer  and  expressed  as  [76], 

Ci1  = - eo<^r(E h)  for  cylindrical  pores  of  radius  R0  14a 


»(r>- -[»,-»(£ -8 .  +  ")!  +  P) 

Differentiation  of  Equation  (7)  with  respect  to  r  yields  the 
following  boundary  condition  at  the  Stern/diffuse  layer  interface  at 
r=R0  +  H, 


ch  _  eoer(En)  ft  for  spherical  pores  of  radius  R0  1 4b 


3.4.  Method  of  solution  and  data  processing 


-e0er{EH)  ^p(r  =  R0  +  H) 


eoer(Eff)  Ro 
RoIog(l+H/R0)Ro+H 

x[is-i(r  =  R0  +  H)} 


Equation  (1)  was  solved  using  the  commercial  finite  element 
solver  COMSOL  4.2,  along  with  the  boundary  conditions  given  by 
Equations  (3c)  and  (10)  or  (13)  and  field-dependent  permittivity 
er(£)  given  by  Equation  (5).  The  simulations  were  run  on  a  Dell 
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workstation  Precision  690  with  two  2.66  GHz  Quad-Core  Intel  Xeon 
CPUs  and  64  GB  of  RAMs. 

Due  to  charge  conservation  [Equation  (1)],  the  total  amount  of 
charges  Q  stored  near  the  electrode  surfaces  As  is  equal  to  that 
present  at  the  Stern/diffuse  layer  interface  denoted  by  Ad.  Then,  it 
can  be  computed  by  integrating  the  surface  charge  density 
(eofrE-n)  along  or  Ad  and  expressed  as  [8,77], 

Q  =  J  e0er(E)E-n  dA  =  J  e0er(£)E-n  d A  (15) 

As  Ad 

where  E  =  -V\f/  is  the  local  electric  field  vector,  n  is  the  local 
outward  normal  unit  vector  at  the  electrode  surface  As  or  at  the 
electrode/electrolyte  interface  Ad.  Then,  the  overall  gravimetric 
capacitance  was  estimated  as  [33], 


where  p  and  V  are  the  density  and  volume  of  the  amorphous  carbon 
electrode,  respectively.  The  diffuse  layer  areal  capacitance  Cf  and 
the  total  areal  capacitance  Cs  were  respectively  estimated  as 
[32,33], 

C?  =  A  =  Jr  f  e0er(£)E  n  dA  and  Cs=7^-  (17) 

s  \fdAd  xfdAd  J  rw  i Ms 

At 

Note  that  the  total  areal  capacitance  Cs  can  be  also  equivalently 
calculated  using  the  one-dimensional  series  formula 
1  /Cs  =  1  /Cf4  + 1  /C,p  when  the  Stern  layer  thickness  is  much  smaller 
than  the  electrode  or  pore  diameter  as  considered  here.  In  addition, 
the  Stem  layer  areal  capacitance  Cf4  was  given  by  the  Helmholtz 
model  [Equation  (11)  or  Equation  (14)]. 

The  numerical  convergence  criterion  was  defined  such  that  the 
maximum  relative  differences  in  the  predicted  capacitances 
Cf4,  Cf ,  Cs,  and  Cg  was  less  than  1.5%  when  decreasing  the  mesh  size 
by  a  factor  two.  The  mesh  size  was  the  smallest  at  the  Stern/diffuse 
layer  interface  due  to  the  large  electric  potential  gradient  and  then 
gradually  increased.  The  maximum  mesh  size  was  specified  to  be 
about  0.1  nm  at  all  Stern/diffuse  layer  interfaces  and  remained  less 
than  2.5  nm  anywhere  else  in  the  computational  domain.  The  total 
number  of  finite  elements  was  on  the  order  of  107  for  the  simula¬ 
tions  of  bimodal  ordered  carbon  CP204-S15  with  1  /12th  of  one  unit 
cell  shown  in  Fig.  2(b). 

Finally,  the  numerical  tool  was  validated  based  on  two  equilib¬ 
rium  cases  reported  in  the  literature.  First,  the  equilibrium  electric 
potential  profile  in  the  diffuse  layer  predicted  by  solving  the  MPB 
model  was  validated  against  the  exact  solution  for  planar  elec¬ 
trodes  [6,8,58]  with  er  =  78.5,  ca  =  0.01  and  0.001  mol  L”1, 
v  =  0,  and  \f/D  =  0.1  V.  Second,  the  computed  capacitances  for  the 
Stem  and  diffuse  layers  obtained  from  the  MPB  model  were  vali¬ 
dated  against  their  theoretical  formula  assuming  constant  elec¬ 
trolyte  permittivity  [34,55,56]  for  fs  =  2  V,  c„  =  1  mol  L  ',  and 
a  =  0.68  or  1.40  nm.  Good  agreement  was  obtained  between  our 
results  and  reported  electric  potential  profiles  [6,8,58]  or  capaci¬ 
tances  [34,55,56]  for  all  cases  considered. 

4.  Results  and  discussions 

4.1.  Validity  of  the  new  boundary  conditions 

The  new  boundary  conditions  [Equations  (10)  and  (13)]  were 
used  to  compute  the  capacitances  for  a  single  cylindrical  electrode 
particle  and  a  cylindrical  pore  by  solving  Equations  (lb),  (3c)  and 
(10)  or  Equation  (13),  respectively.  The  parameters  were  chosen 


such  that  a  =  0.68  nm,  c„  =1.0  mol  L-1,  and  —  2  V  with  field- 
dependent  electrolyte  permittivity  er(£)  given  by  Equation  (5).  The 
results  were  compared  with  those  obtained  by  simulating  both  the 
Stern  and  diffuse  layers  in  the  electrolyte  domain  using  the 
conventional  boundary  conditions  given  by  Equation  (3).  Fig.  3 
shows  the  Stern  layer,  diffuse  layer,  and  total  areal  capacitances 
predicted  using  these  two  approaches  as  a  function  of  particle  or 
pore  radius  R0  ranging  from  2.5  to  60  nm.  The  relative  difference  in 
the  values  of  Cf4,  Cf ,  and  Cs  predicted  using  these  two  approaches 
was  less  than  0.2%  for  all  cases  considered  here.  In  addition,  it  is 
evident  that  the  predicted  Cf4  was  much  smaller  than  Cf  for  all 
cases  considered.  Therefore,  the  double  layer  capacitance  Cs  was 
dominated  by  Cf 4  which  is  consistent  with  the  conclusion  drawn  in 
our  previous  study  [34]  when  accounting  for  field-dependent 
electrolyte  permittivity.  Note  that  for  particle  or  pore  radius 
larger  than  40  nm,  the  predicted  capacitance  reached  a  plateau  of 
Cs  =  10.2  pF  cm-2  corresponding  to  that  of  planar  electrodes  as 
discussed  in  Refs.  [32-34,49,50]. 

Similarly,  the  predicted  capacitances  Cf4,  Cf ,  and  Cs  using  these 
two  approaches  were  compared  for  a  spherical  electrode  and 
a  spherical  pore  with  various  radii  ranging  from  2.5  to  60  nm. 
Excellent  agreement  was  observed  in  all  cases.  Overall,  these 
results  demonstrate  that  the  Stern  layer  can  be  accurately 
accounted  for  by  using  the  new  boundary  conditions  given  by 
Equation  (10)  or  Equation  (13)  for  cylindrical  and  spherical  elec¬ 
trodes  or  pores  without  explicitly  simulating  the  Stem  layer  in  the 
electrolyte  domain.  Note  that  the  total  number  of  finite  elements 
decreased  by  about  30%-60%  when  using  the  new  boundary 
conditions  for  simulating  a  single  cylindrical  and  spherical  elec¬ 
trode  or  pore  with  radius  ranging  from  2.5  to  60  nm.  The  corre¬ 
sponding  computational  time  was  reduced  by  about  10%— 30%.  This 
reduction  in  finite  elements  and  computational  time  became  more 
significant  with  increasing  radius  and  geometric  complexity.  This 
clearly  demonstrates  the  advantage  of  the  new  boundary 
conditions. 

4.2.  Capacitances  of  ordered  bimodal  carbons 

The  double  layer  capacitances  Cf4,  Cf ,  and  Cs  of  the  ordered 
bimodal  carbon  shown  in  Fig.  2(b)  were  predicted  by  solving 
Equation  (lb)  in  the  diffuse  layer  along  with  the  new  boundary 
conditions  given  by  Equations  (3c)  and  (10)  or  Equation  (13).  The 
electrode  surface  was  divided  in  three  sections  (i)  the  inner  surface 
of  the  pores  of  radius  95  nm,  (ii)  the  outer  surface  of  radius  97  nm, 
and  (iii)  the  mesopore  with  7  nm  in  radius  located  in  the  walls 
separating  the  macropores.  The  boundary  condition  given  by 
Equation  (10)  was  imposed  at  the  outer  surfaces  while  Equation 
(13)  was  imposed  at  the  mesopore  and  inner  pore  surfaces, 
respectively.  Note  that  without  these  new  boundary  conditions,  it 
was  impossible  to  solve  the  coupled  governing  equations  for  such 
a  complex  domain  due  to  (i)  the  difficulty  in  creating  the  geometry 
and  (ii)  an  excessively  large  number  of  finite  elements. 

Table  1  summarizes  the  predicted  Stern  layer,  diffuse  layer,  and 
total  areal  capacitances  as  well  as  the  gravimetric  capacitance  at  the 
inner,  outer,  and  mesopore  surfaces  of  the  ordered  bimodal  carbon. 
Results  were  obtained  based  on  the  non-solvated  effective  ion 
diameter  a  =  0.68  nm  [73]  with  field-dependent  electrolyte 
permittivity  er(£)  given  by  Equation  (5).  The  areal  capacitances 
predicted  for  planar  electrodes  for  the  same  parameters  were  also 
reported  in  Table  1  for  comparison  purposes. 

It  is  evident  that  the  Stern  layer  areal  capacitance  Cf4  was  about 
one-third  smaller  than  the  diffuse  layer  areal  capacitance  Cf  at  all 
surfaces.  Thus,  the  total  areal  capacitance  Cs  was  controlled  by  the 
Stern  layer.  In  addition,  the  total  areal  capacitances  Cs  at  the  inner 
and  outer  surfaces  were  10.2  pF  cm-2  which  was  identical  to  that 
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Table  1 

Predicted  Stem  layer  Cfc,  diffuse  layer  C?,  and  total  Cs  areal  capacitances  as  well  as 
gravimetric  capacitance  for  the  ordered  bimodal  carbon  (Fig.  2(b))  computed  at  the 
inner,  outer,  and  mesopore  surfaces.  Results  were  obtained  by  solving  Equations 
(lb),  (3c)  and  (10)  or  Equation  (13)  using  non-solvated  or  solvated  ion  diameter 
a  =  0.68  or  a  -  1.40  nm,  respectively  along  with  c*  -  1  mol  L_1,  \j/s  =  2  V,  and 
electrolyte  permittivity  given  by  Equation  (5).  Predictions  for  planar  electrodes  are 
reported  for  comparison. 


Ion  diameter 

Capacitance 

Planar  Inner 

electrode  surface 

Outer 

Mesopores 

a  =  0.68  nm 

1  Cf'(pFcm-; 

!)  13.4 

13.3 

13.4 

12.8 

C?  (pFcm-2 

)  43.2 

44.6 

43.8 

47.3 

CsCpFcnr2: 

)  10.2 

10.2 

10.2 

10.1 

Cg  (Fg-1) 

— 

39.0 

35.0 

23.0 

a  =  1.40  nm 

l  C|l  (M Fcm-; 

!)  13.0 

11.8 

11.9 

10.0 

C?  (pFcm“2 

:)  17.7 

21.4 

20.7 

28.5 

CsfixFcm-2: 

)  7.5 

7.5 

7.5 

7.4 

Cg  (F  g-1) 

- 

29.2 

25.1 

17.8 

for  planar  electrodes.  Indeed,  electrodes  with  radius  of  curvature 
larger  than  40  nm  behave  like  planar  electrodes  as  established  in 
Refs.  [32,34]  and  also  shown  in  Fig.  3.  Similarly,  the  areal  capaci¬ 
tance  Cs  at  mesopore  surfaces  was  10.1  pF  cm-2  which  fell  within 
1%  of  that  of  planar  electrodes.  In  fact,  it  was  found  that  the  areal 
capacitance  Cs  of  mesopores  decreased  with  increasing  pore  depth 
(or  carbon  wall  thickness)  (not  shown).  This  was  due  to  the 
confinement  of  the  electric  field  in  small  pores  leading  to  reduced 
surface  charge  density  and  capacitance  as  previously  discussed 
[32,33],  However,  this  effect  appeared  to  be  negligible  for  small 
pore  depth  (i.e.,  carbon  wall  thickness)  of  2  nm  such  as  those 
considered  here.  In  addition,  the  mesopores  contributed  about  30% 
less  to  the  gravimetric  capacitance  Cg  than  the  inner  and  outer  pore 
surfaces.  This  was  due  to  the  small  carbon  wall  thickness  and  thus 
the  relatively  small  mesopore  surface  area  compared  with  those  of 
inner  and  outer  pore  surfaces. 

Table  1  also  summarizes  the  predicted  capacitances  for  the 
electrode  structure  shown  in  Fig.  2(b)  but  using  solvated  ion 
diameter  a  =  1.40  nm  [75].  It  demonstrates  that,  for  larger  effec¬ 
tive  ion  diameter,  the  Stern  layer  areal  capacitance  Cf1  had  the  same 
order  of  magnitude  as  the  diffuse  layer  areal  capacitance  Cf  at  all 
surfaces.  Consequently,  the  capacitances  Cf*  and  Cf  contributed 
nearly  equally  to  the  total  areal  capacitance  Cs.  Here  also,  the  inner 
and  outer  pore  surfaces  contributed  slightly  more  to  the  gravi¬ 
metric  capacitance  than  mesopores.  However,  both  the  local  areal 
and  gravimetric  capacitances  decreased  by  about  30%  when  using 
the  solvated  ion  diameter  a  =  1.40  nm  instead  of  a  =  0.68  nm. 
This  was  due  to  the  associated  reduction  in  the  maximum  ion 
concentration  cm  =  1  / NAa 3  at  the  electrode  surface  as  ion  diam¬ 
eter  a  increases. 

Finally,  for  the  relatively  large  surface  potential  fs  =  2  V 
considered,  the  ion  concentration  given  by  Equation  (2)  reached  its 
maximum  value  cm  at  all  surfaces  in  the  simulations  (not  shown) 
regardless  of  the  value  chosen  for  the  effective  ion  diameter  a. 

4.3.  Comparison  with  experimental  data 

Table  2  shows  the  predicted  overall  areal  and  gravimetric 
capacitances  Cs  and  Cg  defined  by  Equations  (16)  and  (17)  and 
accounting  for  the  contribution  of  all  electrode  surfaces.  The 
experimental  values  reported  in  Ref.  [15]  are  also  reproduced  for 
comparison.  It  is  worth  noting  that  the  values  of  Cs  and  Cg  predicted 
using  non-solvated  ion  diameter  a  =  0.68  nm  were  about  8.5%  and 
3%  larger  than  their  respective  experimental  counterparts.  On  the 
other  hand,  Cs  and  Cg  predicted  using  solvated  ion  diameter 
a  =  1.40  nm  underestimated  the  experimental  values  by  about 
24%.  In  fact,  the  predictions  using  non-solvated  and  solvated  ion 


Table  2 

Predicted  overall  areal  and  gravimetric  capacitances  of  the  ordered  bimodal  carbon 
using  (a)  non-solvated  ion  diameter  a  =  0.68  nm  and  (b)  solvated  ion  diameter 
a  =  1.40  nm  in  comparison  with  experimental  data  reported  in  Ref.  [15],  Results 
were  obtained  by  solving  Equations  (3c),  (lb)  and  (10)  or  Equation  (13)  with 
Coo  =  1  mol  L-1  and  =  2  V. 

Specific  area  (m2  g-1)  Cs  (pF  cnr2)  QtFg-1) 
Measured  [15]  1003  04  95 

a  =  0.68  nm  960  10.2  97.9 

a  =  1 .40  nm  960  7.5  72.0 


diameters  set  the  upper  and  lower  bounds  for  the  capacitances, 
respectively.  The  capacitances  predicted  using  non-solvated  ion 
diameter  (a  =  0.68  nm  [73])  showed  better  agreement  with 
experimental  data.  We  speculate  that  these  results  could  be 
attributed  to  two  possible  reasons.  First,  as  the  electrolyte  concen¬ 
tration  increases,  the  dissolved  electrolyte  ions  become 
less  solvated,  i.e.,  they  are  surrounded  by  less  solvent  molecules 
[16,78—80],  Consequently,  the  effective  ion  diameter  decreases  with 
increasing  electrolyte  concentration  [16,78,79],  Note  also  that  the 
solubility  of  TEABF4  in  propylene  carbonate  is  about  1  mol  L  1  at 
room  temperature  [81].  Second,  the  effective  ion  diameter  of  TEA+ 
tends  to  decrease  under  large  local  electric  field  [82],  In  fact,  TEA+ 
ions  was  found  to  become  distorted  and  able  to  adsorb  in  pores  with 
diameter  even  smaller  than  the  non-solvated  TEA+  [82-84],  Thus, 
the  effective  ion  diameter  near  the  electrode  surfaces  approaches 
that  of  the  non-solvated  ions,  i.e.,  a  =  0.68  nm  [82— 84],  Overall,  the 
predicted  capacitances  agreed  well  with  experimental  data.  These 
results  validate  the  numerical  models,  boundary  conditions,  and 
constitutive  relationships  developed  here  for  simulating  EDLCs  with 
three-dimensional  ordered  structures.  The  new  boundary  condi¬ 
tions  were  essential  in  obtaining  such  results. 

Thanks  to  this  experimentally  validated  numerical  model,  it 
becomes  possible  to  numerically  explore  the  effects  of  the  electrode 
architecture  on  its  energy  storage  capabilities.  Fig.  4  shows  the 
predicted  gravimetric  capacitance  Cg  of  bimodal  carbon  structures 
as  a  function  of  their  specific  surface  area  ranging  from  910  to 
1030  m2  g-1.  Here,  the  specific  surface  area  was  varied  by  changing 
the  inner  pore  radius  R0  from  50  to  150  nm  while  other  geometric 
parameters  such  as  the  carbon  wall  thickness  and  mesopore  radius 
remained  identical  to  those  of  CP204-S15  carbon  shown  in  Fig.  2(b). 


Fig.  4.  Predicted  and  experimentally  measured  [15]  gravimetric  capacitance  Cg  for 
bimodal  carbons  as  a  function  of  their  specific  surface  area.  Numerical  results  were 
obtained  by  solving  Equation  (lb)  using  boundary  conditions  given  by  Equations  (3c) 
and  (10)  or  Equation  (13)  with  non-solvated  ion  diameter  a  =  0.68  nm, 
Coo  =  1  mol  IT1,  ^  =  2  V,  and  the  electrolyte  permittivity  given  by  Equation  (5).  The 
inner  pore  radius  R0  was  varied  from  50  to  150  nm. 
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The  results  were  obtained  by  solving  Equation  (lb)  subject  to 
boundary  conditions  given  by  Equations  (3c)  and  (10)  or  Equation 
(13)  using  both  solvated  and  non-solvated  ion  diameters,  i.e., 
a  =  1.40  nm  and  a  =  0.68  nm,  respectively.  Other  parameters 
were  identical  to  those  used  to  generate  the  results  presented  in 
Table  1.  Fig.  4  also  shows  the  measured  gravimetric  capacitance  of 
different  bimodal  carbon  films  [15]  obtained  using  1  mol  L-1 
TEABF4  electrolyte.  It  is  evident  that  predicted  and  experimentally 
measured  gravimetric  capacitance  Cg  increased  linearly  with 
increasing  specific  surface  area.  The  slope  of  Cg  vs.  specific 
surface  area  corresponds  to  a  constant  areal  capacitance  of 
Cs~7.4  or  10.2  fiF  cm-2  when  using  solvated  or  non-solvated 
ion  diameter,  respectively.  Note  that  this  trend  is  consistent  with 
experimental  data  reporting  a  linear  relationship  between  gravi¬ 
metric  capacitance  and  specific  surface  area  with  an  areal  capaci¬ 
tance  of  9.4  pF  cm-2  [15].  Overall,  very  good  agreement  was 
observed  between  experimental  measurements  and  model 
predictions. 

Finally,  it  is  evident  that  the  new  boundary  conditions  [Equa¬ 
tions  (10)  and  (13)]  developed  here  made  possible  the  simulations 
of  EDLCs  with  three-dimensional  ordered  electrode  structures. 
These  simulations  can  give  detailed  information  such  as  the  local 
charge  storage,  electric  potential,  and  ion  concentrations  within  the 
electrolyte  which  cannot  be  measured  experimentally.  Note  also 
that  these  boundary  conditions  [Equations  (10)  and  (13)]  can  be 
readily  employed  to  simulate  the  dynamic  charging  and  discharg¬ 
ing  of  EDLCs  with  ordered  electrode  structures.  Then,  the  model 
could  be  used  to  identify  the  optimum  electrode  architecture  and 
provide  design  rules  to  achieve  maximum  charging  performance  by 
EDLCs. 

5.  Conclusion 

This  paper  developed  a  three-dimensional  (3D)  model  based  on 
continuum  theory  for  simulating  EDLCs  with  ordered  electrode 
structures.  For  the  first  time,  a  new  set  of  boundary  conditions  was 
derived  to  account  for  the  Stern  layer  forming  near  planar,  cylin¬ 
drical,  and  spherical  electrodes  as  well  as  cylindrical  and  spherical 
pores.  They  made  possible  the  simulations  of  EDLCs  with  3D 
ordered  electrode  structures  while  simultaneously  and  accurately 
accounting  for  (i)  both  the  Stern  and  diffuse  layers,  (ii)  finite  ion 
size,  and  (iii)  the  dependency  of  electrolyte  permittivity  on  the  local 
electric  field.  The  model  was  used  to  faithfully  simulate  an  actual 
EDLC  consisting  of  complex  3D  ordered  bimodal  carbons  in 
1  mol  L1  TEABF4/propylene  carbonate  electrolyte  solution  [15]. 
The  predicted  gravimetric  capacitance  of  different  bimodal  carbons 
was  found  to  increase  linearly  with  increasing  specific  surface  area 
corresponding  to  constant  areal  capacitance.  Numerical  predictions 
were  in  very  good  agreement  with  experimental  data  [15]. 
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